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ABSTRACT 

Polytropic models play a very important role in galactic dynamics and in the theory 
of stellar structure and evolution. However, in general, the solution of the Lane-Emden 
equation can not be given analytically but only numerically. In this paper we give a good 
analytic approximate solution of the Lane-Emden equation. This approximation is very 
good for any finite polytropic index n and for the isothermal case at a level < 1%. We 
also give analytic expressions of the mass, pressure, temperature, and potential energy as 
a function of radius. 
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1. Introduction 



As stellar encounters are little important for galaxies, clusters of galaxies, or globular 
clusters (Binney & Tremaine 1987), the fundamental dynamics describing these systems 
is that of a collisionless system. In the general case the collisionless Boltzmann equation 
cannot be solved because it involves too many independent variables. However, we can 
get certain exact solutions of the collisionless Boltzmann equation for a subset of possible 
stellar-dynamical equilibria. The system with a polytropic state equation, corresponding 
to isotropic velocity dispersion tensors, is one of them. In stellar structure and evolution 
theory, the polytropic model also plays an important role (Chandrasekhar 1939). 

For a polytropic system, the relation of pressure P and density g is given by 

P = Kg 1 = Kg 1+ ™ , (1) 

where K is the polytropic constant, 7 is the adiabatic index, and n the polytropic index. 
K is fixed in the degenerate gas sphere (e.g. in a white dwarf or in a neutron star) and free 
in a non-degenerate system. In galactic dynamics n is larger than ^ (Binney & Tremaine 
1987) which means that no polytropic stellar system can be homogeneous, n ranges from 
to 00 in the case of the theory of stellar structure and evolution (Kippenhahn & Weigert 
1989, Chandrasekhar 1939). With the polytropic relation (1), hydrostatic equilibrium, and 
Poisson's equation for gravitational potential field, we derive the Lane-Emden equation 
(see, for instance, Binney & Tremaine 1987, Kippenhahn & Weigert 1989) for spherical 
symmetry: 

d 2 u 2 dw _ 

ie + iTi~- u ■ (2) 

The dimensionless variables u> and £ are defined by 

C = Ar , A 2 = —— — —Q c n , 

(n + 1)K 

N l/n (3) 

.Qc, 

where g c is the density at the center of the sphere and G the Newtonian gravitational 
constant. 00 corresponds to the dimensionless gravitational potential. We have excluded 
the isothermal case n = 00, which we will discuss in detail in a later section. 

The Lane-Emden equation must be solved with the original central conditions: 

KW-i , g) £=o = o , (4, 

which will ensure the regularity of the solution at the center. The solution gives w as a 
function of £. From uj we can get the density profile g. Only for the three values of n = 
0, 1, and 5 can the solution be given in analytic form. Apart from these three cases the 
Lane-Emden equation has to be solved numerically. As the Lane-Emden equation has a 
regular singularity at £ = 0, we can expand a>(£) in a power series as 
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From (2) and (4) we have f^f = -f = 2 ^ at £ = 0. In this paper, we make the 
approximation of taking the second derivative term as to obtain a good approximate 

analytic solution of the Lane-Emden equation. We give our analysis for the isothermal 
Lane-Emden equation in section 2, comparing our results with numerical solution. The 
general case is dealt with in section 3 followed by discussion and conclusion in section 4. 

2. Lane-Emden equation for isothermal sphere 

The Lane-Emden equation for an isothermal sphere (see Kippenhahn & Weigert 1989, 
Binney & Tremaine 1987) can be written as 

where uj is related to the mass density at dimensionless radius £ by 

Q = Qc e~^ , (7) 

and q c is the density at the origin. It is easy to find (see also Binney & Tremaine 1987) 
that 
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w = -ln^J , (8) 

is one of the solutions of Eq. (6). This solution describes a model known as the singular 
isothermal sphere. Unfortunately, the singular isothermal sphere has infinite density at 
£ = 0. To obtain a solution that is well behaved at the origin, equation (6) has to be 
integrated with the central conditions 

U (0) = , (^=0 . (9) 

Its solution can not be given by analytic expressions but only numerical computation. 
The Lane-Emden equation has a regular singularity at £ = 0. In order to understand the 
behaviour of the solutions there, a power series expansion similar to (5) can be derived and 
has to be used. For large radius where the effect of the central conditions is very weak the 
solution should asymptotically approach the singular isothermal solution. The isothermal 
sphere consisting of an ideal gas has an infinite radius as well as an infinite mass. 

At small radius, a useful approximation to g(£) is the modified Hubble law (Binney 
& Tremaine 1987), which was introduced empirically by King (1962), 



g\ 1 



(10) 



Comparing this relation with a numerical solution of (6) one can say that for £ < 5 the 
relative error is less than 5 %. Expression (10) does not fit the isothermal profile well at £ 



3 



> 9 as it approaches asymptotically to a logarithmic slope —3 and not —2, as is required 
by the isothermal profile. 

Since (j^j = 1/3 = | (f^f) at £ = from (6), it is reasonable to approxi- 
mate the second derivative term in (6) with 



1 /2duA 



as our first approximation and to get a first differential equation 



2£d£ e ' UJ 



which can be integrated immediately, giving 



Wl = In ( 1 + |- ) , (12) 



and 

01 



= e~ Wl = w , (13) 



1+S 



where we have used the condition (9) to determine the integration constant. The subscript 
1 indicates a first approximation. 

When £ is very large, Eq. (12) becomes 



Wl ~ In ( L ) , (14) 



which has a very similar profile to (8) and the same tangent as the numerical solution. 
Figure 1 shows the numerical solution of equation (6) and our first approximation (13), 
which shifted by 1/3 approximately reaches the singular solution for large £. 
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Fig. 1 The first order (dashed line indicated by Qi/g c ) and the second order approximation 
(dashed^hne^ indicated by Q2I Qc ) 'of the' density for the isothermal case.' The first orcler 
approximation ^has been multiplied by l/3lePshow its same tangent as thQ-fiumerical result 
(solid line) for large £. £ 



Taking the second derivative of u\ as our approximation to the second derivative term 
and replacing the right hand side term in equation (6) with (13), we obtain the second 
approximation to the solution of equation (6) 

u 2 = 2\n^l + ^-(^-^ + con S , (15) 

or 

— = e 2 = cons 1 H e V / (16) 

Qc V 6 / 

where subscript 2 indicates the second approximation. Integration constant cons in (15) is 
determined by the central condition (9) and is zero. Equations (16) as an approximation 
to the solution of (6) is good for small £. For large £, the approximation has an asymptotic 
behavior, g/g c oc £ -4 , which does not fit the isothermal profile as the same reason for the 
modified Hubble law discussed above. We also show the second approximation in figure 1. 
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As the first order approximation (12) is good for large £ and the second order approxi- 
mation (15) is good for small £, we combine them and obtain a more general approximation 



ujq = auj2 + (1 — a)u>i 

= (1 + a) In (l + -L_^_ + cons . ( 17 ) 

where < a < 1, coming from the condition (9), which requires uj = at £ = 0. For very 
large £, equation (17) asymptotically approaches 

/£ 2 \ 

(1 + a) In I — J — a + cons, 
but not equation (8). So, we change the constant term in (17) to 

to make equation (17) approach (8) for very large f . From (8) and (17) A takes 3 1 / Q 12e. 
Defining the general approximation as u = ujq + Auj and substituting it in equation (6), 
we get the final analytic approximation 



w = (l + a)ln( 1 + C " ° 



6 J 6 i + Si 

1 ^ 6 



n In ( I f " : - 



3 1 /«12e 



) +to (l + (2--l) T ^ I ) 



(18) 



and 

Q 



( 1+ 3T7^) a(^/6)/(l+4) 



(19) 



where e is the natural constant, and D a constant depending on a. Equation (19) reaches 
the singular solution for £ — > oo. The discrepancy of (19) to the numerical solution is 
therefore at its largest value when £ is intermediate. Adjusting a and D, we can modify 
the fitness for intermediate £. One good combination is a = 0.551 and D = 3.84 * 10~ 4 , 
which makes the largest relative error of equation (19) to the numerical solution of (6) be 
0.72%. In figure 2, we compare the structure of densities obtained both numerically and 
analytically and show the relative error (g n — g) / { g n , where g n is the numerical computation 
solution of (6). For £ — > oo, the relative error approaches zero. As in the isothermal case 
u has no physical meaning, we give only the result for density g. The largest error for u 
is also about 0.72%. 
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Isothermal Case 
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3. Lane-Emden equation for general case 

We can apply a similar method for the general case. We first discuss the n^l case. 
The limit case n = 1 will be discussed later. For 1, equation (2) also has a singular 
solution as in the isothermal case 



,2 



(1-nY t2 



l/(l-n) 



"=%rt)«'l ' ,20) 

When n > 5, the solution of equation (2) under central condition (4) should asymptotically 
approach (20) for £ — > oo. We approximately solve Eq. (2) with the same technique as for 
the isothermal sphere. We approximate the second derivative term in Eq. (2) with §|% 
and get 

2 + 52^ _ 

— ' (21) 
where 5 is a constant. 5 can affect the fitness of the approximation strongly for n < 1 
but not for n > 1. We can however get very good approximation if fixing 5 and adjusting 
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other parameters (see below) for any n. Therefore, we will always take 5 = 1. Integrating 
equation (21), we get 

u 1 = (l + A n e) 1/il - n) , (22) 

where A n is (n — 1) /6. Here, we have used the central condition (4) to determine the 
integration constant. However, this first analytic approximation does not have the right 
behaviour for large £ compared to the numerical solution of equation (2). Equations (22) 
for any n > 1 can not reach zero at finite £ and give an infinite radius. Integrating Eq. (2) 
gives 

«=-/W-/i$« ■ ^ 

Substituting (22) in the right-hand-side of Eq. (23) and integrating, we get our second 
approximate analytic solution: 

U2 = cons + 2 (1 + A n i 2 ) ^ + ^ (1 + A n i 2 ) ^ , (24) 

where the constant cons is determined by the central condition (4) as —1. This equation, 
however, reaches zero at finite £ for any finite n and gives a radius smaller than the one 
obtained from numerical computation. 

Since solution (22) gives too large values and (24) gives too small values for large £, 
we construct a more general approximation for w asa linear combination of equations (22) 
and (24) 

ujq = auj2 + (1 — a)u>i 

= -a+(l + a) {1 + A n t 2 ) ! + (1 + A n f) 

We would like to point out that if we keep 5 in (21) as a free parameter, then equation 
(25) gives a good approximation to the solution of (2) at level < 1% with proper 5 and a 
for n < 1. In this case 5 varies from to 1. We notice that (25) asymptotically reaches 
the constant —a not zero when £ — > oo for n > 5. As (25) should approaches the singular 
solution at very large £ (see also, Eggleton, 1995) as in the isothermal case, we change the 

constant —a in Eq. (25) to —a (l + B^ 2 ) 1 ^ 1 " and let (25) for £ — > oo equal to the 
singular solution (20) when n » 5 and to zero when n = 5. We get 

BiH ^"-'(i+^-(§^r~T" • ; (26) 

n(n-l) 6 ( 4a \ r , r 

where we have considered the conditions that B n should equal zero at n = and 1 and 
continue at n = 5 to get the expression for n < 5. With B n given by (26), (25) gives a 
good approximation with relative error less than 5% for n < 5 and about 1% for n » 5. 



8 



n=1 .5 



0.6 



0.4 



Solid: Numerical result 
Dashed: Analytical approximation 



CO 



0.2 



Fig. 3 The analytic approximation before a Auj modification given by equation(20)4s^faixly 
good. Here is the result for n = 3/2 as an example. For n = 3/2 case, the relative erroris 



less than abcsut 1% when a is taken as 0.44 
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Figure 3 gives the results for n = 1.5 as an example. For n = 1.5, the relative error is less 
than about 1% when a = 0.44. 
So, if our general solution is 



to = too + Ato , 



(27) 



Ato would be very small. From Eqs.(2) and (25) as well as (26), we get at very large £ for 
n = 5 and at £ — > £ n for n < 5, 

Aw~- , (28) 
as where u>o ~ 0. In order to avoid the singularity of Auj at £ = 0, we set 



Ato -— 

(1 + D n ^f 

where C n and D n are small constants and the index (3 determined from Eq.(2) is 
=6.47 - 7.01/3i + 5.53/3^ - 25.63/3 2 + 49.42/^ - 26.88/3| , 

9 



(29) 



(30) 
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where ft = l/(l + (n- 5) 2 ) and (3 2 = l/(l + (n- 3) 2 ) . With the same procedure and 
the same reason as for n < 5, we have for n > 5 



1 + F n e? (i + Dn £(jy 



A ^ ~ i , F ^ + ~ — ^,2 > ( 31 ) 



where the index rj is (n — 2) / (n — 4) and the coefficient E n determined from equation (2) 
is 



5 — n / n — 4 7i — 4' 



with which our approximation always asymptotically approaches the singular solution (20) 
at large £ for n > 5. To fit the solution for intermediate £, the parameter F n can be given 

as 

f (™-5) 6 f ^ e 

F n ~ <( (n-3) 5 (4n+50)' tOT " " 5 5 (33) 

0, for n < 5 . 



Finally, we configure a general result for any n 



= - a (1 + B n £ 2 ) 



2\V(l-n) 



+ (l + a) {l + A n e) 



2\l/(l-»») 



+ f£ 2 (1+A n £ 2 ) 



(34) 



where 77 is (n — 2) / (n — 4), A n is (n — 1) /6, B n is given by (26), E n by (32) and (3 by 
(30). The free parameters a, F n , C n and D n should be determined by fitting the numerical 
solution of equation (2). (34) with proper values of the parameters a, C n , D n and F n can 
give a very good approximation solution of equation (2) under the central condition (4). 

For n » 5, equation (34) tends to the singular solution as £ — > 00. In this case, (34) 
always approximates the solution of Eq.(2) with < 0.1%. Especially when n > 10, the last 
term and when n > 500 the last two terms in equation (34) are not important and can 
be ignored and the approximation is within the level 0.1%. As n decreases, (34) slowly 
deviates the singular solution to finally reach uj = 1/(1 +£ t 2 /3) 1 / 2 at £ very large for n = 5 
where the sum of the first three terms in (34) approach zero in accordance with (25) and 
(26). (34) reaches zero at infinity for n > 5 and at finite £ n for n < 5. All the properties of 
(34) discussed above are also the properties of the exact solution of (2) (a detail dicussion 
to the general properties of the solution of equation (2), see Binney & Tremaine (1987), 
Kippenhahn & Weigert (1989)). 
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(34) is the exact solution for n = 0. For general n, we have to adjust a, C n , D n 
and F n to get a good approximation. Equation (34) gives a good approximation to the 
solution of equation (2) for small and large £, which is insensitive to the parameters within 
proper ranges of values. For intermediate £, however, special a, C n , D n and F n are needed 
to reduce the error of the approximation. In figures 4 and 5, we show the variations of 
potential uj, density p and relative error {uj n — uj)/u n with radius £ for n = 3/2 and for 
n = 5, respectively. For n = 3/2, the relative error is larger near the finite radius £ n and 
oscillates around the mean value zero. Such oscillation behavior also exists in figure 5 for 
n = 5. In fact, the oscillation always exists for any n. We give the values of the parameters 
used to get figures 4 and 5 and the largest relative fitting errors e max in table 1. In table 1, 
we give the values of the parameters for several polytropic indices n and the largest fitting 
errors e max . We also give the best fitting value for F n in the table. With the values of 
a, C n and D n and the value of F n given by equation (33), equation (34) can also give an 
approximation at level e max < !%• 
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For n < 5, the relative error at £ — > £ n is larger as the numerical solution u n approaches 
zero. (34) with B n given by (26) reaches zero at finite £ and is the exact solution for n = 0. 
For n > 5, the error is at its maximum when £ is in the range of around 1 to 20, depending 
on index n. When £ is very large, the approximation approaches the singular solution. 
With the increasing distance to the center of the gas sphere, the effect of the central 
boundary (4) on the structure of the sphere decreases. Therefore, the solution under the 
central boundary condition (4) approaches the singular solution. The approximation is 
very good even for £ — > oo. For n > 10, C n can be set to zero. For n > 500, E n is less 
than 10 -3 and Au given by equation (31) becomes unimportant. 

For n = 1, integrating (21) and considering the central boundary condition (4) we get 
the first approximation 

wi = e~ ei% . (35) 
Then, with the same procedure as for general n, we get our second approximation 
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Table 1. Parameters for several polytropic indices n 

n polytropic index, a mixing parameter of the first and the second approximations, e max 
(given in percent) maximum relative error of analytic approximation to numerical calculation 
result. For other parameters, see the text. 



n 


a 








IIL LI JL \ 


0.5 


0.5 


1.30659 * 10 4 


2.3 * 10~ 4 





0.37 


1 


0.455 


1.27746 * 10 4 








0.15 


1.5 


0.481 


4.61841 * 10 4 


11 i n— S 

1.1 * 10 13 





0.12 


2 


0.512 


o on 1 o i r\ — 4 

8.3218 * 10 


2.56 * 10 





0.28 


3 


0.53 


5.56215 * 10" 4 


2.745 * 10" 2 





0.60 


5 


0.545 


1.9415 * 10" 3 


3.348 * 10" 2 





0.34 


6 


0.56 


2.0* 10" 4 


3.0* 10" 2 


5.3* 10" 5 


0.50 


7 


0.526 


1.53* 10" 7 


l.l * 10" 3 


7.40* 10" 4 


0.39 


8 


0.540 


6.1 * 10" 8 


1.0* 10" 3 


2.7* 10" 3 


0.35 


9 


0.546 


3* 10" 8 


9* 10" 4 


5.4* 10" 3 


0.26 


10 


0.538 


1.12 * 10" 8 


9* 10" 4 


8.91 * 10" 3 


0.24 


20 


0.613 







4.2 * 10" 2 


0.12 


50 


0.65 







o.ioi 


0.060 


500 


0.58 







0.19 


0.021 



= _l + 2e -« 2 /6 + i! e -C 2 /6 . ( 36) 



Wo 
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(35) and (36) are also the limit of expressions (22) and (24) for n — * l, so the uj is continuous 
at n = l. From the limit case of equations (26) and (34) for n = l, we therefore get the 
final result for the n = l 

u) = -ae to I t+wz ) « 4. (i + a ) e — _|_ -Ce~ H , , (37) 

6" (l + D^/ 3 ) 2 

where /3 is given by equation (30), C n and D n are constant. In table l, we also give a set 
values of C n and D n and the largest approximate error of equation (37). For n = l, the 
largest relative error of the approximation is only 0.15%. 

4. Discussion and conclusion 

We gave our analytic approximation solutions to the Lane-Emden equations for isothcr-| 
mal sphere in section 2 and for general polytropic cases in section 3. From the analytic 
results we can derive some useful relations. 

From (3) and (2) for the general case and from (7) and (6) for the isothermal case, 
the mass contained in a sphere of radius r 



m = I Aixr gdr 

(38) 



o 



., , 1 du 

4w 1 ^ 
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where — sign is for the general case and + for the isothermal case. Dimensional r and 
non-dimensional £ are related for the general case by equation (3) and for the isothermal 

case by 

For the general case and n^l, it follows from (34) that 



m(r) =4rcg c r 3 

(^(HiU 1 )*^!!^)*^^^ , (4Q) 
_ f „ +T% _ 2 1+T^ + T^^" _ „ rfa _ 3 2^-l-g„^ 

For n < 5, the star has a finite radius R or £ n , which is estimated from (34) by setting 
cj to zero, and a finite mass M. From (40) we can get the expression for mass M from 
which the central density g c is determined. However, for n > 5 the polytropic gas sphere 
has infinite radius and infinite mass. In this case we cannot estimate the central density 
directly. From equations (l)and (3) we have 

n+l 

P = Kg c n iv n+1 , (41) 

where to is given by (34) and K is determined by equation (3) as A is given by A = ^ n /R 
for n < 5. For an ideal gas the temperature is given as 



kg k 



Kg l J n u , (42) 



where \i is the mean molecular weight, m v is the mass of atomic Hydrogen, and k is the 
Boltzmann constant. From the definition of potential energy within radius r of the sphere 

E g = - ^-dm> , (43) 
Jo r 

the polytropic relation (1), and the hydrostatic equilibrium equation, we get for n ^ 5 

3 Gm 2 3 n+l. 3n .... 

E q = $m - Anr 6 P , (44) 

5 — n r 5 — n 5 — n 

where $ is gravitational potential, given by 

§ = -(n + l)Kg^ . (45) 

This means that for n < 5 the polytropic sphere has finite radius and the potential is set 
to zero at the surface. For n > 5 the potential becomes zero at infinity. 
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In polytropic models, n = 3 and n = 3/2 are two important cases. We have given n = 
3/2 as one example before. Here, as another example, we construct a polytropic model of 
index 3 of the sun (M = 1.989* 10 33 g and R = 6.96* 10 10 cm) and compare our results with 
that from numerical computation. For n = 3, from table 1 and equations (26) and (30), we 
have a = 0.53, (3 = 2.199, B n = 4.648 * 10" 3 , C n = 5.56215 * 10" 4 and D n = 2.745 * 10" 2 . 
Then we get £ 3 = 6.897 from (34) and A = £ 3 /R = 9.895 * 10 -11 . Substituting these in 
equation (40) we get central density g c = 76.96 g cm" " 3 , K = 3.87* 10 14 from equation (3) 
and consequently, P c = 1.26* 10 17 dyn/cm 2 from equation (41) as £ = at the center. For 
the ideal gas \i = 0.62 we get from equation (42) the central temperature T c = 1.2 * 10 7 K. 
From Kippenhahn and Weigert (1990), numerical computation gives £3 = 6.897, g c = 76.39 
g cm -3 , P c = 1.24 * 10 17 dyn/cm 2 and T c = 1.2 * 10 7 K. We see that our approximation 
for n = 3 gives a quite good result at a level 0.7% for density, 1.6% for pressure and 0% 
for temperature. 

For the isothermal sphere, the mass from equations (18) and (38) is 
m =8irg c r s 

fl + a 6a a D (2~ a — 1) \ (46) 

1^6 + £ 2 " (6 + £2)2 " 3iA>i 2e + ^ + (1 + 2~ a D^ 2 ) (1 + D^ 2 ) J ' 

where a = 0.551 and D = 3.84 * 10" 4 . When £ » ^6 the mass changes with radius 
as m oc r. The pressure is given by polytropic relation (1) for n = 00. The constant 
K is related to the central density g c by (39). However, as the isothermal sphere does 
not have a finite radius, the constant A cannot be determined as in the case n < 5. In 
the isothermal sphere the temperature is constant everywhere. From definition (43) and 
the hydrostatic equilibrium equation, we find the potential energy within radius r for the 
isothermal sphere 

E g = 4nr 3 Kg - 3Km , (47) 

where g is given by (19) and m by (46). 

In this paper we have given a good analytic approximate solution of the Lane-Emden 
equation with which we have obtained analytic expressions for the mass contained in radius 
r, the pressure, temperature, and gravitational potential energy within radius r. It would 
be interesting to apply the approximation in modeling the structures of stars and galaxies. 
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